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ABSTRACT 

We calculate the one-point probability density distribution functions (PDF) and the power 
spectra of the thermal and kinetic Sunyaev-Zeldovich (tSZ and kSZ) effects and the mean 
Compton Y parameter using the Magneticum Pathfinder simulations, state-of-the-art cosmo¬ 
logical hydrodynamical simulations of a large cosmological volume of (896 Mpc//i) 3 . These 
simulations follow in detail the thermal and chemical evolution of the intracluster medium as 
well as the evolution of super-massive black holes and their associated feedback processes. 
We construct full-sky maps of tSZ and kSZ from the light-cones out to z = 0.17, and one 
realisation of 8°.8 x 8°.8 deep light-cone out to z — 5.2. The local universe at z < 0.027 
is simulated by a constrained realisation. The tail of the one-point PDF of tSZ from the deep 
light-cone follows a power-law shape with an index of —3.2. Once convolved with the effec¬ 
tive beam of Planck, it agrees with the PDF measured by Planck. The predicted tSZ power 
spectrum agrees with that of the Planck data at all multipoles up to l tv 1000, once the 
calculations are scaled to the Planck 2015 cosmological parameters with ii rn = 0.308 and 
as = 0.8149. Consistent with the results in the literature, however, we continue to find the 
tSZ power spectrum at l = 3000 that is significantly larger than that estimated from the high- 
resolution ground-based data. The simulation predicts the mean fluctuating Compton Y value 
ofY = 1.18 x 10" 6 for n m = 0.272 and cr 8 = 0.809. Nearly half (« 5 x 1CT 7 ) of the 
signal comes from halos below a virial mass of 10 13 Mr.fiIt. Scaling this to the Planck 2015 
parameters, we find Y = 1.57 x 10~ 6 . 

Key words: hydrodynamics, method: numerical, galaxies: cluster: general, cosmic back¬ 
ground radiation, cosmology: theory 


1 INTRODUCTION 

One of the last unsolved problems in cosmology with the cosmic 
microwave background (CMB) is connected with unknown levels 
of the spectral distortion of the black-body spectrum of CMB. Two 
decades ago the Far Infrared Absolute Spectrophotometer (FIRAS) 
on board NAS A’s Co smic B ackgro und Explorer (COBE) provided 
the first limits l IFixsen et al.lll996ll on the so-called y- and /7-type 
spectral distortions, which originate from energy releases in th e 
present and early universe ISunvaev & Zeldovicblll970bl . Il980al l. 
An enormous progress in technology of cryogenics and detec¬ 
tors of millimetre and sub-millimetre radiation enables us now to 
reach sensitivity which is 100 or even 1000 times better than that 
of FIRAS. A rece ntly proposed PIXIE instrument is an example 
dKogut et al.ll2014h . 


* E-mail: dolag@usm.uni-muenchen.de 


While measuring the absolute monopole of the spec¬ 
tral distortions requires an absolute spectrometer such as 
FIRAS and PIXIE, the fluctuating component does not. 
Kh atri & S unvae v (2015) have recently dem onstrated that the 
Planck data i Planck Collaboration et al.ll2015d) permit to limit the 
mean of the fluctuating part of the y-type distortion to be Y < 
2.2 x 10~ 6 (where Y is the so-called Compton Y parameter), 
thereby reducing the observational upper bound by a factor of seven 
compared to the original COBE/FIRAS limitQ 

The y-type spectral distortions do not carry information about 
redshifts at which they were produced. Therefore, it will be a chal¬ 
lenge to distinguish between the primordial y-distortions originated 
from the pre-recombination era and those from a low-redshift uni- 


1 However, we should not forget that the COBE/FIRAS limit applies to 
the sum of the uniform and the fluctuating parts, whereas the Planck limit 
applies only to the latter. 


© 0000 RAS 














2 K. Dolag, E. Komatsu, R. Sunyaev 


verse induced by the large-scale structures. Nevertheless, in the 
standard thermal history of the universe, the latter component is 
expected to dominate by a couple of orders of magnitude]; namely, 
the dominant component of the mean Y is expected to come from 
the fluctuating component of the thermal Sunyaev-Zeldovich (tSZ) 
effect, i.e., inverse Compton scattering of low energy CMB photons 
on hot non-relativistic electro ns in groups and clusters of galaxies 
dSunvaev & Zeldovicblll970ah . Therefore, the new Planck bound 
provides a limit on the total thermal energy content of hot gas in 
the universe. Indeed, this bound is approaching the predicted value 
of the mean Y calculated by the previous generation s of cosmo ¬ 
logical hydrodynamical simulations; for example, iRefregier et al.l 
feOOCh find Y = 1.67 x 1CP 6 for a flat A Cold Dark Matter (CDM) 
model with = 0.37 and ag = 0.8 (also see lda Silva et alj200d ; 
[Seliak et aljEotlllzhang eTaI1l2004bi) . It is therefore timely to re¬ 
visit this calculation with much improved, state-of-the-art simula¬ 
tions. 

The Planck satellite has measured the power spectrum of fluc¬ 
tuations of tSZ on large angular scales (l < 1000). The tSZ power 
spectru m is very sensitive to the amplitude of matter density fl uctu- 
ations JlComatsu & Kitavam31l999l ; lKomatsu & Seliak|l2002h . The 
power spectrum on large angular scales is particularly a powerful 
probe of the amplitude of fluctuations, as it is less sensitive to as- 
trophysics of the core of galaxy clusters than tha t at l « 3000 
l l Komatsu & Kitavamall 19991 : iMcCarthv et alj|2014h . 

Cosmological hydrodynamical simulations have proven es¬ 
sential in interpreting the observational data on the statistics 
of tSZ such as the power spectrum and the one-point proba¬ 
bility density distribution function (PDF). The simulations have 
steadily improved over the past two decades, in terms of the 
numerical methods, the volume, the mass and spatial resolu¬ 
tion, and implementation of baryonic physics suc h as cool¬ 
ing, heating, chemistry, and feedback processes i Persi et all 


19951: IRefreeier et all200d:lda Silva et al.l200Cl:ISeliak et all200ll; 


SjrringdetaL 2001 1 _ daSilvaetaL 2001b|: lRefreg ier& Te^ssierl 

' lWhite et al.l 120021; IZhang et al.l 120021; iDolag et al.l l2005tJ: 

; jRoMarefli^talJJ^OOTljshaw^etal w, 

120121 : Munshi et aiT 201 3 I : McCarthy et al.1 


Hallman et al.l 

2007) 

Battaglia et al. 

2010 


The recent progress in performing larger and more “complete” 
cosmological hydrodynamical simulations, which follow baryons 
in different p hases such as in Warm-Hot Intergalactic Medium 
(WHIM lCen & Ostriketl dl999l) ). Intracluster Medium (ICM), stars 
and super massive black holes and their related feedback, allows 
now for a direct and detailed comparison of the predicted tSZ sig¬ 
nals with observations. Such comparison may shed light on the 
recently reported tensions between cosmological parameters in- 
ferred_from the pri mary C MB with the ones inferred from tSZ 
dPlanck Collaboration et al l2015d lbh. 

Scattering of CMB photons off electrons moving with a non¬ 
zero line-of-sight bulk velocity with respect to the frame of co- 


2 In principle, knowing the precise expectation for the low-redshift tSZ 
contribution would make it possible to estimate the detectable level of the 
primordial y- type distortions from the pre-recombination era. In the ab¬ 
sence of primordial non-Gaussianity, we do not expect any angular depen¬ 
dence of the y- type signal before recombination. Therefore, characterizing 
the angular dependence of the fluctuating component of the y-signal from 
the large-scale structures might open a possibility to separate that contribu¬ 
tion from the primordial monopole (via, e.g., cross-correlation of tSZ with 
other tracers of the large-scale structure). If successful, we can retrieve in¬ 
formation on the energy release at redshifts of 1000 < z < 20000. 


ordinates where the CMB is isotropic yields temperature fluc¬ 
tuations by the Doppler shift of light, and this _is known as 
the kinetic Su nyaev-Zeldovich (kSZ) effect dSunvaev & Zeldovichl 
Il970ain~980bh . There are two contributions to the kSZ signal: the 
reionisation contribution from 2 > 6, and the post-reionisation 
contribution. The calculation of the former requires detailed reion¬ 
isation simulations including radiative transfer and is beyond the 
scope of this paper. The calculation of the latter is in principle sim¬ 
pler than the former. Precisely characterising the post-reionisation 
kSZ is important, as we must subtract this contribution from the 
total kSZ to extract the reionisation contributi on w hi ch, i n turn, 
can constrain th e physics of reionisation (see IZahn eta] .1120121 ; 
Park et aT] |2()l 3l . and references therein). Roughly speaking, the 
post-reionisation kSZ power spectrum is twice as large as that from 
reionisation; thus, ten per cent uncertainty in the post-reionisation 
contribution resultsjn twenty per cent uncertainty in the reionisa¬ 
tion contribution I Park et all2015t) . 

In a fully ionised universe, the kSZ effect is given by the line- 
of-sight integration of the radial momentum field, i.e., a large-scale 
velocity field mod ulated by a small-scale density fluctuations of 
electrons (see lPark et a] .II20T5I . and references therein). Therefore, 
the calculation of the post-reionisation kSZ requires a simulation 
whose box size is large enough to capture the long-wavelength 
velocity field, and the spatial resolution high enough to resolve 
non-linear structures of baryons responsible for scattering of CMB 
photons. The hydrodynamical simulations of the post-reio nisation 
kSZ have also improved over the past decade dSnringel et all 


KaZ nave also improved o ver the past decade ( spnngelet_alj 

200 ll; | da Silva et al .11 2001 alibi: I White et all|2002l: IZhang et al.ll200i 


2004ts; Roncarelli et al.ll2007l : Shaw et afll20 1 2l : Dolag & Sunvaevl 

20131) . The measur ements of the post-reionisation k SZ are improv¬ 
ing rapidly as well l lHand et all2012l : ILi et al.ll20l4h . 

In this paper, we shall push the simulation frontier on the in¬ 
vestigation of tSZ and kSZ further. The paper is organized as fol¬ 
lows. Section 2 reviews the simulation method and the light-cone 
generation. Section 3 shows detailed comparisons of the simula¬ 
tion results with the observational data on the Coma cluster, and the 
one-point PDF and power spectrum of tSZ. We also show the sim¬ 
ulation predictions for kSZ toward Coma, and the one-point PDF 
and power spectrum of kSZ. In section 4, we study how the mean 
Compton Y signal builds up over cosmic time. We summarise our 
findings in section 5. 


2 SIMULATIONS 


We construct sky maps of tSZ and kSZ using two sets of simu¬ 
lations: the “local universe” simulation (Sec. I2H based on a con¬ 
strained realisation of the local universe at 2 < 0.027, and the Mag- 
neticum Pathfinder simulation (Sec. l2.2j . Combining these simula¬ 
tions, we construct full-sky maps of tSZ and kSZ out to 2 = 0.17, 
and one realisation of 8°.8 x 8°.8 deep light-cone out to 2 = 5.2. 

Both simulations are based on the parallel cosmological Tree 
Particle-Mesh (PM) Smoothed-particle Hydrodynamics (SPH) 
code P-GADGET3 dSpringell l2005h. The code uses an entr opy- 
conserving formulation of SPH ( Springel & Hemquistl 12002 ) and 
follows the gas using a low-visc osity SPH s cheme to properly tr ack 
turbulence foolag et alj|2005bl) . Based on iDolag et all d2004l) . it 
also follows thermal conduction at l/20th of the classical Spitzer 
value dSnitzerl 19621) . It also allows a treatment of radiative cooling, 
heating from a uniform time-dependent ultraviolet background, and 
star formation with the associated feedback processes. 

Radiative cooling rates are computed by following the same 
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procedure presented by [Wiersma et al.l I 2 OO 9 I) . We account for 
the presence of the CMB and the ultraviolet (UV)/X-ray back¬ 
ground radiatio n from quasars and galaxies, as computed by 
lHaardt & Madaul ( i200ll) . The contributions to cooling from each 
one of 11 elements (H, He, C, N, O, Ne, Mg, Si, S, Ca, Fe) have 
been pre-comput ed using the publicly available CLOUDY photo¬ 
ionization code dFerland et al .11 199~ij) for an optically thin gas in 
(photo-)ionization equilibrium. 

We model the interstellar medium (ISM) by usin g a s ub¬ 
resolu tion model for the multiphase ISM of ISpringel & Hernauistl 
( 2003;). In this model, the ISM is treated as a two-phase medium, in 
which clouds of cold gas form by cooling of hot gas, and are em¬ 
bedded in the hot gas phase assuming pressure equilibrium when¬ 
ever gas particles are above a given threshold density. The hot gas 
within the multiphase model is heated by supernovae and can evap¬ 
orate the cold clouds. A certain fraction of massive stars (10 per 
cent) is assumed to explode as supernovae type II (SNII). The re¬ 
leased energy by SNII (10 51 erg) triggers galactic winds with a 
mass loading rate proportional to the star formation rate (SFR) with 
a resulting wind velocity of tVnd = 350 km/s. 

We include a det ailed model of chemical evolution accord¬ 
ing to iTornatore et ah . ( 2007 ). Metals are produced by SNII, by 
supernovae type la (SNIa) and by intermediate and low-mass 
stars in the asymptotic giant branch (AGB). Metals and energy 
are released by stars of different masses, by properly account¬ 
ing for m ass-dependent life-times (w ith a lifetime function ac¬ 
cording to lPadova ni & Matteuccill 1 9931). the metallicity-dependent 
ste llar yields by Woos lev & Weavei 119951) for SNII, the yields 
by Ivan den H oek & Groeneweeenl 1 1 9971) for AGB stars, and the 
yields by Thielemann et al. 1 20031) for SNIa. Stars of different 
masses are init ially distr ibuted according to a Chabrier initial mass 
function (IMF: lChabrieill2003l) . 

Most importantly, our simulations include prescriptions for 
the growth of black holes and the feedback fro m active galac¬ 
tic nuclei (AGN) based on the model of ISpringel et alj d2005l) 
and Di^atteojstaJ 1 20051) with the same modifications as in 


iFabian et al.ld201()l) ~and some new, minor changes as described be 
low. 

The accretion onto black holes and the associated feedback 
adopts a sub-resolution model. Black holes are represented by col¬ 
lisionless “sink particles”, which can grow in mass by either ac¬ 
creting gas from their environments, or merging with other black 
holes. The gas accretion r ate, M ,, is estimated by the Bondi-Hoyle- 
Lyttleton approxim ation dHovle & Lvttletor]| 1 9391 : [Bondi & Hovlel 
1 19441: iBondill 1952h : 


M. = 


4ttG 2 M? /boostp 
(c§ + l; 2)3/2 ’ 


(1) 


where p and c B are the density and the sound speed of the sur¬ 
rounding (ISM) gas, respectively, /boost is a boost factor for the 
density which typically is set to 100 and v is the velocity of the 
black hole relative to the surrounding gas. The black hole accretion 
is always limited to the Eddington rate, i.e., the maximum possible 
accretion for balance between inwards-directed gravitational force 
and outwards-directed radiation pressure: M, = min(M., M e dd). 
Since the detailed accretion flows onto the black holes are unre¬ 
solved, we can only capture black hole growth due to the larger 
scale, resolved gas distribution. 

Once the accretion rate is computed for each BH particle, the 
mass continuously grows. To model the loss of this gas from the 
gas particles, a stochastic criterion is used to select the s urround- 
ing gas particles that are accreted. Unlike in lSpringel et al.l ( 120051) . 


in which a selected gas particle contributes to accretion with all its 
mass, we include the possibility for a gas particle to accrete only 
with a fraction (1/4) of its original mass. In this way, each gas par¬ 
ticle can contribute with up to four ‘generations’ of BH accretion 
events, thus providing a more continuous description of the accre¬ 
tion process. 

As for the feedback, the radiated luminosity, L r , is related to 
the black hole accretion rate by 

L r = e r M.c 2 , (2) 

where e r is the radiative efficiency, for which we adopt a fixed 
value of 0.1. This value is assumed typically for a radiatively 
efficient accretion disk onto a non-rapid l y spinning black hol e 
dShakura & Sunvaevll 19731 : ISpringelll2005l : |Pi Matteo et alJl2005l) . 
We assume that a fraction ef of the radiated energy is thermally cou¬ 
pled to the surrounding gas; thus, E{ = e r efM.c 2 is the rate of the 
energy feedback, ef is a free parameter and typically set to 0.1. The 
energy is distributed to the surrounding gas particles with weights 
similar to SPH. In ad dition , we incorporate the feedback prescrip¬ 
tion of IFabian et al.l d2010t) : namely, we account fo r a transition 
from a quasar- to a radio-mode feedback (see also ISiiacki et al J 
I2OO7I) whenever the accretion rate falls below an Eddington-ratio 
°f fe dd = M./M e dd < 1(T 2 . During the radio-mode feedback 
we assume a 4 times larger feedback efficiency than in the quasar 
mode. This way, we attempt to account for massive black holes that 
are radiatively inefficient (having low accretion rates) but are effi¬ 
cient in heating the ICM by inflating hot bubbl es in corresp ondence 
of the termination of AGN jets. In contrast to lSpringel et al .1 d2005l) . 
we modify the mass growth of the BH by taking into account the 
feedback, e.g., AM. = (1 — t r )M.At. We introduced some more 
technical modifications of the^ original implementation, for which 
readers can find details in lHirschmann et al] d2014f) . where we also 
demonstrate that the bulk properties of the AGN population within 
the simulation are quite similar to the observed AGN properties. 


2.1 Local Universe Simulation 


The local universe simulation uses the final output of a cosmolog¬ 
ical hydrodynamical simulation of a constrained realisation of the 
local uni verse based on the nearly full-sky IRAS 1.2-Jy galaxy sur¬ 
vey data jFisher et al.ll 1994 , 1995|). T he ini tial conditions are simi¬ 
lar to those adopted bv lMathis et al] (1 20021) in their study of struc¬ 
ture formation in the local universe. The galaxy distribution in the 
IRAS 1.2-Jy galaxy survey is first smoothed by a Gaussian with 
a width of 7 Mpc and then linear ly evolved back in time up to 
z = 50 following the method of iKolatt et alj d 19961) . The result¬ 
ing fie ld is th en use d as a constraint on phases of random Gaussian 
fields dHoffman & Ribakil 199 ll) for initialising the simulation. 

The volume constrained by the observational data covers 
a sphere of radius ~ 80 Mpc /h centered on the Milky Way. 
This region is sampled with more than 50 million high-resolution 
dark matter particles and is embedded in a periodic box of ~ 
240 Mpc//i on a side. The region outside the constrained volume 
is filled with nearly 7 million low-resolution dark matter particles, 
allowing a good coverage of long-range gravitational tidal forces. 
The gravitational force resolution (i.e., the softening length) of the 
simulation has been fixed to be 7 kpc /h (Plummer-equivalent), 
fixed in physical units from z = 0 to z = 2 and then stays constant 
in the corresponding co-moving units (e.g. 21 kpc /h) at higher 
redshifts. _ 

Unlike in the original simulation made bv lMathis et alJd2002h . 
where only the dark matter component is present, here we follow 
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Figure 1. Full-sky maps of the Compton Y parameter (Eq.[3| top) and the kSZ effect (Eq.[4| middle) obtained from the local universe simulation (z < 0.027) 
combined with the full sky maps of the Magneticum Pathfinder simulation (0.027 < £ < 0.17). We also show 8° .8 x 8° .8 maps from the deep light-cone of 
the Magneticum Pathfinder simulation restricted to 0.17 < z < 5.2. The lower panels show, from the left to right: the Compton Y from the deep light-cone 
(0 < £ < 5.2) smoothed with a 9.66 arcmin FWHM Gaussian beam; the Compton Y at the native resolution of HEALPix with nside=2048; the kSZ at the 
native resolution, and the sum of the two at 150 GHz (Eq.[5). 
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Figure 2. Compton Y from the deep light-cone (0 < z < 5.2). The right panel shows a zoom onto a region containing several rich clusters at various redshifts. 


also the gas and stellar components. Therefore, we extend the ini¬ 
tial conditions by splitting the original high-resolution dark mat¬ 
ter particles into gas and dark matter particles having masses of 
m Bas « 0.48 x 10 9 Mq/A and mj m ~ 3.1 x 10 9 tf@/A, respec¬ 
tively; this corresponds to a cosmological baryon fraction of 13 per 
cent. The total number of particles within the simulation is slightly 
more than 108 million, and the most massive cluster is resolved by 
almost one million particles. 

In this simulation, we assume a flat ACDM model with a 
present matter density parameter of = 0.3, a Hubble constant 
of Ho = 100 h km/s/Mpc with h = 0.7, and an r.m.s. density 
fluctuation of as = 0.9. 


2.2 Magneticum Pathfinder Simulation 

The Magneticum Pathfinde simulation follows a large 
(896 Mpc/A) 3 box simulated u sin g 2 x 1526 3 particles and 
adapting a WMAP7 j Komatsu et al .1120 171) flat ACDM cosmology 
with erg = 0.809, A = 0.704, Sl m = 0.272, f4 = 0.0456, and an 
initial slope for the power spectrum of n s = 0.963. Dark matter 
particles have a mass of rriDM = 1-3 x 10 10 M@/h, gas particles 
have m gas ~ 2.6 x 10 9 M@/ft depending on their enrichment 
history, and stellar particles have m sta , ra ~ 7.5 x 10 s JV/q/A 
depending on the state of the underlying stellar population. For 
gas and dark matter, the gravitational softening length is set to 
10 kpc/A (Plummer-equivalent), fixed in physical units from 
2 = 0 to 2 = 2 and then stays constant in the corresponding 
co-moving units (e.g. 30 kpc/A) at higher redshifts. For star 
particles it is accordingly half the values (e.g. 5 kpc/A at 2 = 0). 
In this simulation, one gas particle can form up to four stellar 
particles. 


www. magneticum .org 


2.3 Map Making 

We construct maps from the simulations using SMAC l lDolag et al .1 
l2005al) . integrating both the tSZ and kSZ signals through our hy- 
drodynamical simulations. 

The tSZ signal in each pixel at a sky position 9 is characterised 
by the Compton Y parameter defined by dSunvaev & Zeldovichi 

IT9721) 

Ytsz(d) = [ d l n B (0, l) T(d, l) , (3) 

m e c 2 J 


where n e and T are the three-dimensional number density and tem¬ 
perature of thermal electrons, respectively, and Ab, ctt, m e , and c 
are the Boltzmann constant, the Thomson scattering cross section, 
the electron mass, and the speed of light, respect i vely. T he kSZ sig¬ 
nal is obtained by dSunvaev & Zeldovlchlll970al . ll980bl) 

w ksz (9) = = -Zr fdl ne (9, l) v r (9, l ) , (4) 

J- cmb C J 

where T C mb = 2.725 K and v r is the radial component of the 
peculiar bulk velocity, defined such that v r > 0 for gas moving 
away from us. 

As the tSZ and kSZ have different dependence on the observed 
frequencies, we use the following weighted sum of the two when 
showing the combined signal: 

vfrq _ g{x)YtS7. — WkSZ , c ^ 

^eff — / \ 5 W/ 

g(x) 

with 


g{x) 


x{e x + 1 ) 

(e x - 1) 


( 6 ) 


derived by IZeldovich & Sunvaevl dl969l ). where x = 
frq[GiT2]/56.8. The temperature anisotropy due to tSZ is 
given by <5Ttsz/Tcmb = g(x)Y, where we ignore relativistic 
corrections. 
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Figure 3. Compton Y profile of the Coma cluster (left panel). The symbols show the Planck data taken from 1 Planck Collaboration et alj i2013cb . These data 
points have an estimate of the mean background (horizontal dashed line) subtracted by the Planck team. The blue line shows the Compton Y toward the 
counterpart of the Coma cluster in the local universe simulation, while the thick black solid line shows the Coma cluster embedded in the deep light-cone 
of the Magneticum Pathfinder simulation, to show the effect of the contribution of gas outside the Coma cluster. All simulated maps are smoothed with a 
Gaussian beam of 10 arcmin. The error bars on the thick black solid line show the 25% and 75% percentile of the Compton Y distribution in each radial bin, 
averaged over 9 different maps obtained by placing Coma at different positions within the light-cone. The horizontal dotted line shows the mean Compton Y 
of the deep light-cone, V = 1.18 X 10 l ' . The magenta line shows the sum of tSZ and kSZ (Eq.|5j at 150 GHz. In the local universe simulation. Coma is 
moving away from us, thus giving a negative kSZ near the center. On the other hand, a positive kSZ in the outskirt is due to a massive, gas rich sub-structure 
moving toward us, as shown in the mass-weighted velocity map of Coma in the local universe simulation (right panel). 


We produced_ full-sky mags of the Compton Y and kSZ in 
the HEALPix dGorski et alJI 19981) format with Nside=2048 from 
the local universe (0 < z < 0.027) and from the Magneticum 
Pathfinder simulation covering 0.027 < z < 0.17. We also pro¬ 
duced one realization of a deep, 8°.8 x 8°.8 light-cone covering 
0 < z < 5.2. When we combine the deep light-cone with the full- 
sky maps, we use only the relevant parts, 0.027 < z < 5.2 or 
0.17 < z < 5.2, of the deep light-cone. 

In figure |T| we show the full-sky maps of the Compton Y 
(top panel) and kSZ (middle) where the local universe simulation 
is combined with the Magneticum Pathfinder simulation. We also 
show 0.17 < z < 5.2 of the deep light-cone (as indicated by the ar¬ 
row) to give the impression of the additional SZ signals covered by 
the deep light-cone. The lower panels show the Compton Y from 
the deep light-cone (0 < z < 5.2) smoothed with a 9.66 arcmin 
FWHM Gaussian beam as well as the native resolution and the 
same with the kSZ signal added at 150 GHz (Eq.[5)- The deep light- 
cone contains about 8000 galaxy clusters and groups with virial 
masses above 10 13 ' 5 M@/h, contributing to the most prominent 
structures visible in the zoom onto a sub-part of the Compton Y 
map shown in the right panel of figure[2] 


3 COMPARISON WITH THE PLANCK DATA 
3.1 Coma 

The thermal electron pressure determines the (non-relativistic) tSZ 
effect. Azimuthally-averaged radial profiles of thermal pressure in 


galaxy clusters identified in our hydrodyna mical sim ul ations fol - 
low the so-called “universe pressure profile” dAmaud et al.l2010al ). 
The stacked pressure profiles of galaxy-cluster-size halos in our 
simulation are in good agreement with the stacked pressure pro¬ 
files infe rred from th e tSZ data on galax y clusters detected by 
Planck JPIanck C ollaboration et al. |2013ai) and South Pole Tele¬ 
scope (SPT) ^McDonald et al Jl2014h . However, such comparisons 
allow only a statistical comparison of the averaged profiles. 


Local, well-resolved galaxy clusters enable a more detailed 
object-by-object comparison. Here we make use of the fact that 
the constrained simulation of the local universe allows to cross- 
identify objects in the simulations with the real-world counterparts. 
Despite the relatively low spacial resolution (e.g., > 5 M pc ) of 
the co nstraints used for initialising the simulation (see lMathis et al.l 
120021 for details), local galaxy clusters like the Coma cluster have 
a remarkably similar counterpart in the simulation. 

In the left panel of figure [3] we compare the radial profiles 
of tSZ (the Compton Y parameter given in Eq. [3}, as well as the 
significant contribution of kSZ at 150 GHz (as given in Eqs.[4]and 
toward the Coma cluste r in the simulation w ith the Planck tSZ 
data for Coma dPlanck Collaboration et al.ll2013d ). The blue line 
shows the Compton Y toward the “Coma cluster” in the local uni¬ 
verse simulation out to a co-moving distance of 110 Mpc, while 
the thick black solid line shows the “Coma cluster” embedded in 
the deep light-cone out to z — 5.2. The latter has more signal in 
the outskirt of the cluster because of the contributions from hot gas 
outside Coma. 


We find an excellent agreement between the simulation and 
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Figure 4. PDFs of the Compton Y parameter, normalised such that f dy P(y) = 1 with y = 10 6 y. (Left) The black line shows the PDF of the deep 
8°.8 X 8°. 8 light-cone out to z = 5.2 without any smoothing applied, while the blue dashed-triple-dotted line shows the PDF smoothed with 1.1 arcmin 
FWHM beam to simulate the signal that would be measured by SPT. The red dashed line shows the PDF including the local universe simulation, smoothed 
with 10 arcmin FWFIM beam to simulate the signal that would be measured by Planck. (Right) Comparison to the Planck data. The black symbols show 
two estimates (“LIL" and “MILCA/NILC”) of PDFs from the Planck data iKhatri & Sunvaevll2015l : lKhatrill2015l) . The “LIL” method yields larger noise in 
the Compton Y map, and thus has a broader core in the middle, while both methods agree well in the tail. The red and blue lines show the PDFs from the 
simulation with Gaussian noise of cr = 1.5 X 10“ 6 (red) and 2.5 X 10 — 6 (blue), which approximate the estimates of noise in LIL and MILCA/NILC maps, 
respectively. The excess at large y values above the Pl anck data is due to th e str uctures like Perseus-Pisces that are masked when the measurement is done. 
Applying the same mask with f s ^ y = 0.51 as used bv lKhatri & Sunvaevl feoist) . we find that the PDFs from the simulation (solid lines) agree well with the 
Planck data at all values of y. 


the Planck data (shown as the symbols) in the inner part, r < 
100 kpc. The full light-cone integration (thick black line) agrees 
reasonably well with the Planck data up to a few Mpc. Given the 
uncertainties in the constrained simulation, this level of agreement 
is remarkable. Also the mean Compton Y (Y = 1.18 x 10~ 6 ) 
obtained from the deep light-cone (shown by the horizontal dotted 
line) agrees with the background subtracted from the Planck data 
(shown by the horizontal dashed line) to within a factor of two. 

Finally, to get a feel for the magnitude of kSZ toward the 
“Coma cluster”, we show the sum of tSZ and kSZ at 150 GHz 
(Eq.HJ in the magenta line. Both include the full light-cone inte¬ 
gration out to 2 = 5.2. For this particular realisation, we find a 
non-negligible (10 to 20 per cent) contribution from kSzQln par¬ 
ticular, we find that the local universe simulation predicts that the 
overall halo of Coma is moving away from us at « 400 km/s with 
respect to the CMB rest frame, yielding a negative kSZ signal up to 
Mpc radius. Being a merging system, the core of Coma in our real¬ 
isation moves with even hi gher velocity (up to 800 km/s) as com¬ 
monl y seen in simulations JZuHone et al.ll20fol : iDolag & Sunvaevl 
120131) near the center. 

We also find a large, positive kSZ signal in the outskirt, which 
is due to a massive, gas-rich infalling sub-structure moving toward 


4 However, we should not compare the magenta line with the Planck data 
in the left panel of figure [3] as the Planck data shown here are obtained by 
combining the multi-frequency Planck data specifically to extract tSZ, and 
most (if not all) of kSZ has been removed. 


us. See the right panel of figure [3]for a map of the mass-weighted 
velocity field around Coma. While the overall halo velocity pre¬ 
dicted for Coma should be accurate to the extent of the precision 
of the constraints used by the local universe simulation, details of 
sub-structures present in this realization of the local universe are 
far outside the predictive power of such a simulation. Therefore, 
the prediction for a positive kSZ in the outskirt and the phase of 
movement of the core should be interpreted with caution. 

3.2 The Compton Y Map 

3.2.1 Simulation results 

In figure [4] we compare the one-point PDF of the Compton Y map 
from the simulations (lines) with that observed by Planck (symbols; 
Planck Collaboration et al.ll2014bl : IKhatri & Sunvaevll2015l : iKhatril 
20151) . The asymmetry of the left and right tails (skewness) of P(y) 
is due to galaxy clusters_and groups along the li nes of s ig ht, as pre- 
dicted previously ( Yoshid a et al .1 1200II : iRubino-M artm & Sunvaev 
!2003i) . In the absence of noise, the pixel values in the simulations 
are always positive, whereas Planck’s Y values can have negative 
values due to noise (or unaccounted other sources). 

In the right panel, we have smoothed the simulated Compton 
Y map with a Gaussian beam with 10 arcmin FWHM, and added 
Gaussian noise to the simulations with a = 1.5 x 10 -6 and a = 
2.5 x 10~ 6 to approximate the noise estimates from “LIL” and 
“MI LCA/NILC” maps of the Compton Y, respectively (see figure 
3 in IKhatri & Sunvaevl 1201 5t IKhatril l2015h . Since Planck cannot 
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measure the absolute value of Y, there is ambiguity in the exact 
zero point of the PDF. While the simulation agrees well with the 
“LIL” map, there seems a slight offset in the zero point with respect 
to the “MILCA/NILC” map. 

As our maps have significantly higher spatial resolution than 
Planck, we find a larger excess of high Y values (black line in the 
left panel) if not smoothed. This is driven by the central parts of 
halos along the lines of sight. This can be best seen in figure [2] es¬ 
pecially in the zoom-in in the right panel. However, once smoothed 
with the beam size of Planck’s Y map (10 arcmin FWHM), the 
excess due to these structures unresolved by Planck is reduced (red 
dashed line in the left panel), and the simulation and the Planck data 
are in much better agreement up to Y ~ 2 x 10 . The effect of 

the beam smearing in the light-cone map can be seen visually in the 
left most panel in the bottom panels of figureQ] With the smoothing 
size as large as this, the excess PDF at larger values of Y above the 
Planck data is dominated by the nearby structures that subtend large 
angles in the sky. We find that the excess is dominated by the struc¬ 
tures in the local universe simulation. Remarkably, we could iden¬ 
tify the source of the excess PDF as the structures such as Perseus- 
Pisces that are masked when the measurement is don e. Usin g the 
same mask that retains 51% of the sky used bv lKhatri & Sunvaevl 
J2015h . we find an excellent agreement between the PDFs from the 
simulation (solid lines in the right panel) and the Planck data at all 
values of Y. The “core” of the PDF from the simulation at small 
values of Y is dominated by the structures in the deep light-cone 
beyond the local universe simulation. 

We also convolve our map with the SPT beam (1.1 arcmin 
FWHM; magenta line in the left panel). In this case, the excess 
at large Compton Y values is only mildly suppressed, demon¬ 
strating that SPT-like instruments would be able to resolve almost 
all the contributions of structures resolved by our simulation. See 
iHill et all J2014h for the measurement and interpretation of the one- 
point PDF of the ACT data. 

Coming back to the unsmoothed PDF, we find that the tail 
of the PDF follows a power law shape over at least two orders of 
magnitude in Compton Y values. The slope of this power law is 
approximately —3.2 as shown by the dashed line in the left panel of 
figure[5] To quantify non-Gaussianity of the PDF, we calculate low- 
order cumulants of y = y x 10 6 . We find 1.184, 4.292, 24.27, and 
1514 for the mean, variance, skewness, and kurtosis, respectively. 

The middle panel shows the PDF of kSZ (solid line) as well 
as the contribution of the 2 = 0 slice. It s hows a non- G aussian 
tail in agreement with the previous w ork ( Ida Silva et al. 1 12001 at 
lYoshida et al.l200ll ; lHallman et alj200^1 . The low-order cumulants 
of w = w x 10 -6 are 0.17, 5.1, 0.054, and 0.71 for the mean, 
variance, skewness, and kurtosis, respectively. The ensemble aver¬ 
age of the mean should vanish; a non-zero value from the simula¬ 
tion is due to cosmic variance. The right panel shows the PDF of 
the sum of tSZ and kSZ (Y e ff defined in Eq. [5} at various Planck 
frequencies. The contribution to the PDF of Y e f f of the kSZ sig¬ 
nal is negligible for Y e g > 10 -5 , while it significantly modifies 
the PDF at smaller values of Y e g. As the PDF of kSZ is flatter at 
small values of kSZ than that of the Compton Y at small values of 
Y, the PDF of the sum of the two shows an excess probability in 
10 ~ 6 1$ Teff Hr 5 (see the solid lines figure[5]i. The dashed lines 
show the PDF for negative values of Y e g. This unique shape of the 
PDF, which changes as a function of frequencies in a predictable 
way, may be used to detect the kSZ signals in the data. 


3.2.2 Analytical model of the PDF of Compton Y 


The tail of the PDF of Compton Y can be calculated analytically 
following lHill et all d2014l) . Specifically, we compute 


r5 d y rM ma . x 

P{Y) = AY dz to dlnM 

2X1 JO dZ JM min 


dn(M , z) 
d In M 


\8 2 (M, z, Y) — 9 2 (M, z,Y + AY)], (7) 


where AY = 10~ 6 , M is the virial mass with (M m i n , M max ) = 
(5 x 10 14 M 0 /h, 5 x 10 15 Mg/ft), dV/dz is the differential 
comoving volume per steradian, and dn/d\nM is the halo mass 
function given by 


dn(M,z) _ dlnM 2 oo m dn(M 2 oom, z) 


dlnM dlnM dlnM 200m ’ 

with dn/dlnM 2 oom derived from simulations. Here, AT 2 oom is 
the mass enclosed within r 2 oo m , in which the mean overdensity 
is 200 times the mean mass density of the universe. We c onver t 
M to M 2 oom using an NFW profilejNavamMrtjiL 19961 fl997h 
with the c oncentration parameter o f iDuffv et al.l 1200 si) (see, e.g., 
Eq. 14 of iKomatsu & Seliakll200ll) . For the mass range we con¬ 
sider, we find d In M 2 oom/d In M f» 1 to good approximation. For 
dn/dlriM 2 oom, we use the mass fu nctio n from the Magneticum 
Pathfinder simulation (Eq. 1 and 3 of lBocquet et al.ll2015l . with the 
parameters for “M 2 oom Hydro” given in their Table 2). 

To compute the linear r.m.s. mass density fluctuation nec¬ 
essary in the _fitting formulae of the mass function, we use the 
CAMB code dkewis et al ]200d) to generate the linear matter power 
spectrum with the cosmological parameters of the Magneticum 
Pathfinder simulation. 

The “0” in Eq.[7]is the inverse function of 


a r 6 i-5oo j - 

Y(9,M 500 ,z) = —E. dlP e (JP + D 2 A 9 2 ), (9) 

m e C J_6r 5 oo V 

where P e (r) is the electron pressure profile and Da is the proper 
angular diameter distance. That is to say, for a given value of Y, 
M 500 , and 2, we find the corresponding value of 9 using this equa¬ 
tion. M500 is the mass enclosed inside rsoo, within which the mean 
overdensity is 500 times the critical density of the universe. Again, 
we convert M to A /5 00 using arpNFW profile with the concentra¬ 
tion parameter of lDuffv et all d2008l) . 

For P e we use the foll owing parametrized profile dNagai et al .1 
l2007l;lArnaud et al .11201 Obi) ; 


P e (x) = 1.65 (h/0.7) 2 eV cm 3 


AF500 


E s/3 (z) 


2/3 -\-otp 


p(x), ( 10 ) 


3 x 10 14 (0.7 /h)M e 
p = 0 . 12 , E(z) = H(z)/H 0 = 


with x = r/rsoo, a. 

[fi m (l + z) 3 + 1 — ll m ] 1// “, and the function p(x) is defined by 

6.41(0.7//i) 3/2 


p(x) = 


(csoo*)t'[ 1 + (c 5 ooo;) a p-T')/ a ’ 


( 11 ) 


with C 500 = 1.81, a. = 1.33, fj = 4.13, and 7 = 0.31 
dPlanck Collaboration et al.ll2013bl) . However, this fitting function 
for the pressure profile was derived for AT500 assuming hydro¬ 
static equilibrium, which is known to be biased low rela tive to the 
true AT 500 due to non-thermal pressure (see, e.g.. Ishi & Komatsul 
1201-4 . and references therein). We thus rescale AT500 as M500 —> 
AT 5 oo/ 1-2 and rsoo as rsoo —► rsoo/l- 2 1//3 to account for the hy¬ 
drostatic mass bias. This correction brings the pressure profiles in 
the Magneticum Pathfinder simulation into good agreement with 
the Planck data dPlanck Collaboration et al.l2013bl) . 
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Figure 5. PDFs of the SZ effects computed from the deep light-cone. (Left) PDF of the Compton Y parameter, which is the same as the black line in figure|4] 
normalised such that f dy P(y) = 1 with y = 10 6 y. The dashed line shows a power-law with a slope of —3.2. The analytical models with M m i n = 2 x 10 13 
and 2 x 10 14 M ^Jh are shown for co mparison. (Middle) PDF of kSZ, normalised such that f dw P(w) = 1 with w = 10 6 iu. The dashed line shows the 
PDF at z = 0 (see lYoshida et all200ll for comparison), while the solid line shows the PDF from the light-cone. (Right) PDF of the sum of the two (Eq.[5j at 
various Planck frequencies, normalised such that f dY e ff P(Y e ff) = 1 with Y e ff = 10 6 Y^ff. The solid lines show PDFs for Y e ff > 0 while the dashed lines 
show those for Y e ff < 0. 


Eq. [7] is valid when halos included in the calculation do not 
overlap i n the sky . The overlapping fraction, F c i ust , as defined by 
Eq. 10 of lHill et al] J2014h . is 0.32 for M min = 2 x 10 14 M @ /h. 
We show the result of the analytical model in the left panel of 
figure [5] The shape of the PDF in the tail is reproduced rea¬ 
sonably well. We also show the result for a smaller mass limit, 
Mmin = 2 x 10 13 Mq/ h, showing contributions from lower mass 
halos. However, this case gives T’ c i ust = 5.5, i.e., haloes overlap 
significantly in the sky, and thus the calculation cannot be trusted 
for lower values of Y. Nonetheless this calculation gives an idea 
about the contributions from lower mass haloes. Since the PDF in 
the extreme tail is sensitive to the core of the pressure profiles, a dif¬ 
ference in details of the pressure profile in the simulation and that 
used in the analytical model shows up there: the simulation gives 
more pressure in the cores of galaxy clusters than that of Eq.llOl 

3.3 The angular power spectra of tSZ and kSZ 

3.3.1 Simulation results 

To calculate the angular power spectrum of tSZ, we separately 
use the full-sky map obtained from the local universe simulation 
in z < 0.027 and the light-cones over 8°.8 x 8°.8 in 0.027 < 
z < 5.2. As these simulations are performed with different cos¬ 
mological parameters, we rescale the amplitudes of the tSZ power 
spectra by erf 12^ to Planck 2015’s best fitting CMB cosmological 
val ues, Q m = 0.308 an d erg = 0.8149 (“TT+lowP+lensing” of 
Ipianck Collaboration et all2015al) . We have checked that this scal¬ 
ing agrees with the scaling predicted by the analytical calculation 
presented in section U.3.21 

Figure[6]shows the tSZ power spectrum measured by Planck at 
150 GHz (red symbols with error bars: |Planck Collaboration et alj 
l2015ch . Two dashed lines show the tSZ power spectrum from the 
local universe at low multipoles and that from the light-cone at high 
multipoles, while the black solid line shows the sum of the two. 
This division in multipoles is consistent with the previous work 
showing that the nearby structures dominate at low multipoles sim¬ 
ply because they appear larger in the sky Refregier et al.l l2QOOl ; 
[Komatsu & Seliakl2002l ; Dolag et a l.l2005al : Hansen et al ■ 20051). 


The black solid line agrees with the Planck data well at 
all multipoles measured by Planck, i.e., I < 1000. Our con¬ 
clusion that the predicted tSZ power spectrum with Planck 
2015’s TT+lowP+lensing parameters agrees with the measured 
power spectrum is consistent with the finding of the Planck 
team JPlanck Collaboration et al.l l2015ch . [McCarthy et alj ]2014l) 
show that the Pla nck 2013 parameters with = 0.3175 and 
as = 0.834 l iPlanck Collaboration et al.ll2014all over-predict the 
tSZ power spectrum, and show that anoth er set of p arameters with 

= 0.302 and as = 0.817 dSpergel et al.1120IT) gives the tSZ 
power spectrum in agreement with the measurement. The latter set 
is close to the Planck 2015 parameters with lensing that we use in 
this paper; thus, our conclusion is consistent with their results. Us¬ 
ing the fi^ag scaling, for example, the former set gives 32% larger 
power than the 2015 parameters, while the latter gives 4% smaller 
power. 

However, at l ~ 3000, our prediction is significantly 

higher than the measurements reported by SPT and Atacama Cos¬ 
mology Telescope (ACT) collaborations l lReichardt et al .1 120121 : 
ISievers et al.l 120131) . This finding is not new: the previous calcu¬ 
lations of the tSZ power spectrum also over-predict the power at 
l ~ 3000 compared to the SPT and ACT data, although the de¬ 
gree of overestimation varies depending on the details of baryonic 
physics implemen t ed in the models (e.g., iMcCarthv et al.l 120141 : 

iGeorge et aT1l20 1 Si: iRamos-Ceia et all 20141) . 

Whether this discrepancy at l ~ 3000 poses a serious chal¬ 
lenge to theory is unclear, given that the SPT and ACT do not have 
as many frequency channels as Planck. Distinguishing the primary 
CMB, tSZ, and the other extra-galactic sources by their frequency 
dependence is more challenging for SPT and ACT. (On the other 
hand, Planck does not have angular resolution to resolve the power 
at l ~ 3000.) The tSZ signal is a sub-dominant contribution to 
the power spectrum at l ~ 3000 compared to the extra-galactic 
sources. The magenta symbols with error bars in figure [6] show 
the SPT power spectrum with the best-fitting primary CMB power 
spectrum subtracted. The difference between the magenta symbols 
and a magenta vertical line at l = 3000 indicates the amount of the 
extra-galactic power that needs to be subtracted. At least, our tSZ 
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power spectrum does not overshoot the magenta data points that 
provide a firm upper limit on the tSZ power at these angular scales. 

Next, we show the kSZ power spectrum from the light-cones 
by the blue solid line in figure [6] Our_result at high multipoles 
(l > 1000) agrees with that of Tshaw et all ( l2012f) . The upturn 
at lower multipoles can be understood by the contribution from 
the longitudinal velocity contributi ons th at did not fully_cancel by 
the line-of-sight integration dl lernandez-Monteagudo & Holl2009l) . 
At 150 GHz, the kSZ amplitude becomes comparable to tSZ at 
l i=3 300, and becomes even dominant at lower multipoles. (Note 
that the Planck data shown in this figure remove most of the kSZ 
signals by construction, and thus should not be compared with the 
blue line.) However, as the kSZ signal is dominated by the largest 
modes present in the simulation (see figureQJ, even larger cosmo¬ 
logical volumes will be needed to obtain a fully converged kSZ 
results from such light-cones. 


3.3.2 Analytical model of the tSZ power spectrum 


To check accuracy of scaling the tSZ power spectrum to other cos¬ 
mological parameters, we compute the tSZ power spectrum using 
an analytical model. Ignoring a small contrib ution from the correla¬ 
tion b etween two distinct dark matter halos dKomatsu & Kitavama 


1999 ), we model the SZ power spectrum as dKomatsu & Seliak 


20021) 


Ci = g (x) 


r dz d j- r 

Jo dz J M 


max .. _ .dn(M. z) N ,2 

d\nM :E \yi(M,z)\ 2 , 


d\nM 


( 12 ) 


where (M m ; n , M max ) = (5x 10 11 M Q /h , 5x 10 15 M Q /h). The 
2D Fourier transform of the Compton Y parameter, yi , is given by 


Vi 


47rr5oo &t 


I 2 

l 500 


m e c‘ 


■ f dxx 2 P e (x) 
Jo 


sm(lx/l 500 ) 

Ix/lsoo 


(13) 


where x = r/r 500 and 1 500 = Da/t 500 - The electron pressure 
profile, P e (x), is given by Eg. 1101 

For dn/dln M 200111 , we use three different sets of fit param¬ 
eters obtained from numerical simulations in the literature. First, 
we use the mass fun ction from the Mag neticum Pathfinder simu¬ 
lation (Eq. 1 and 3 of lBocquet et al.ll2015l . with the parameters for 
“M 200111 Hydro” given in their Table 2. Using “JV^oom DMonly” 
gives a similar result: the difference in the power spectrum is less 
than four percent at all multipoles). This mass function should give 
the result that is most consistent with the black solid line shown in 
figure |6l T he other mass func tion fits are taken from iTinker et al.1 
d2008h and lTinker et all J2010h . 

To compute the linear r.m.s. mass density fluctuation nec¬ 
essary in th e fitting formulae of the mass function, we use the 
CAMB code dLewis et al.ll200ol) to generate the linear matter power 
spectrum with the Planck 2015 “TT+lowP+lensing” parameters: 
fl b h 2 = 0.02226, n c h 2 = 0.1186, Q. v h 2 = 0.00064, h = 
0.6781 , (0.05 Mpc" 1 ) = 2.1 39 x 10' 9 , and n s = 0.9677 

dPlanck Collaboration et al.l2015ah . 

The light blue solid line shows in Fig. [6] the an alytica l model 
calculation with the mass function of lBocquet et al.1 J2015h . which 
is in good agreement with the results obtained directly from the 
simulation at l < 1000, while it under-predicts the power at higher 
multipoles. This can be understood partially as due to inhomogene¬ 
ity in the distribution of pressure within halos, as th e analytical 
model assumes a smooth distribution of pressure. iBattaglia et al.1 
d2012ll show that inhomogeneity increases the power at these high 
multipoles by 20%. The remaining difference may be due to the 



Figure 6. Power spectra of temperature anisotropies due to tSZ and 
kSZ at 150 GHz in units of //.K 2 . The red symbols with error bars 
show the estimatio n o f the tSZ power spectrum from the Planck data 
jpianck Collaborati o n et al"]|2015d) . while the magenta sy mbols with error 
bars show the power spectrum of the SPT data l lReichardt et al.ll20l 2|) with 
the best-fitting primary CMB power spectrum subtracted. The green and 
magenta vertical lines show the ranges o f the estimated tSZ power spectra 
at l = 3000 by ACT Isicvers et al.ll2013h and SPT, respectively. The black 
dashed lines show the tSZ power spectrum of the full-sky, local universe 
simulation (z < 0.027) at low multipoles, and that of the deep, 8° .8 X 8° .8 
light-cone from the Magneticum Pathfinder simulations (0.027 < z < 5.2) 
at high multipoles. The black solid line shows the sum of the two. Both 
power spectra are scaled to Planck 2015’s “TT+lowP+lensing” cosmolog¬ 
ical parameters of fl m = 0.308 and rrg = 0.8149. The light blue lines 
show the analytical predictions based on three different mass functions (see 
section [7.1.21 . The tripple dot-dotted blue line shows the kSZ power spec¬ 
trum of the deep light-cone. 


pressure profile of low-mass halos (less than 10 14 Mq) in the Mag¬ 
neticum Simulation being slightly different from the Planck pres¬ 
sure profile given in Eq. [TO] as well as to the normalisation of the 
mass function for low-mass halos (see the next paragraph). In any 
case, the overall agreement between the simulation and the analyt¬ 
ical model is satisfactory. 

The analytical models with the other mass functions yield 
similar, but differ ent, results. The o nly d ifference between the 
mass functions of ITinker et al.1 d2008h and ITinker et al.1 d20ld) is 
that the latter forces the normalisation of the mass function by 
/ 0 °° d717200m dn/din M 200111 = p, where p is the mean mass den¬ 
sity of the universe. This normalisation mainly changes the abun¬ 
dance of low-mass halos which are not well resolved by their N- 
body simulations. As a result, the latter mass function (shown as 
the dashed light blue line in figure |6) gives larger power at high 
multipoles where the contributions from low-mass halos dominate. 

At l < 100, the tSZ power spectra with both Tinker et al. mass 
f uncti ons are slightl y larger than that with the mass function of 
iBocquet et al.1 d2015l) (20 and 15 per cent larger at l = 10 and 100, 
respective ly). Indeed, the fitting formula for the mass function of 
IBocquet et al .1 d2015l) gives similarly smaller dn/d in Tkfeoom than 
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Figure 7. How the mean Compton Y builds up over time, dY7dln(l + z). 
All redshifts below z K 1,5 contribute to the total signal (black solid line) 
almost equally, to within 30%. The coloured lines show d Y /dln(l + z) with 
high-mass halos (whose virial masses are indicated by the numbers in units 
of Mg / h) removed from the simulation. Half of the signal at z = 0 comes 
from clusters with M > 10 14 Mg /h, whereas at z ft; 1 the bulk of the 
signal comes from lower-mass halos. The dashed lines show the scatter due 
to kSZ at 150GHz. The bulk of kSZ comes from high redshift and non- 
collapsed regions. The gray solid line show approximate estimates of tSZ 
from the epoch of reionisation as w ell as from the inte rgalac tic medium 
(IGM) at lower redshift taken from iKhatri & Sunvacvl <12012 1 . The gray 
dashed line shows the analytical model for the contributions from halos. 


the fits of Tinker et al. mass functions at the relevant mass scales, 
i.e., Afijoom > 10 lj M^. iBocquet et all 1 20151) explain this by the 
way the fits are performed; namely, the actual data of the N-body 
simulations are similar between Bocquet et al. and Tinker et al., 
but the fitting procedures give slightly different results. Bocquet et 
al. use the Bayesian likelihood approach taking into account prop¬ 
erly the Poisson nature of the mass function measured from the 
simulation, whereas Tinker et al. use the \ 2 statistics. The former 
results do not depend on the bin size of the mass in which the fits to 
dn/d In M 2 oom are performed, whereas the latter results do. Our 
results highlight the importance of better understanding the high- 
mass end of the mass function for the study of the tSZ power spec¬ 
trum. 


decreasing redshift like the volume swept by a light-cone of a given 
area in the sky, we produce maps of the full simulation at each red¬ 
shift. This allows us to compute dP/ dz at every snap shot with the 
same precision given by the comoving box size of our simulation, 
which is slightly more than 2 Gpc 3 . 

Figure[7]shows d!7dln(l + z) as a function of z, i.e., the time 
evolution of the contribution per logarithmic redshift interval to the 
overall signal. We find that all redshifts below a ft; 1.5 contribute 
almost equally to within 30%. To study which masses contribute, 
we also show dP/dln(l + z) with high-mass halos above a cer¬ 
tain mass threshold removed from the simulation. At 2 = 0, half 
of the signal comes from M > 10 14 Mq /h, whereas at z ft; 1 
smaller halos dominate. We find that, even at z = 0, there is more 
than 10 % of dP/dln(l + z) coming from outside of resolved ob¬ 
jects, e.g., the diffuse baryon component. This fraction increases at 
higher redshifts and reaches almost 30% at 2 ft; 1. 

We also show an order - of-ma gnitude comparison (gray solid 
line) by IKhatri & Sunvaevi d2012l) of the contribution from the 
epoch of reionisation and from low redshift WHIM. The first one is 
based on the reionisation optical depth inferred from CMB obser- 
v ations b y WMAP whereas the later one is based on the simulations 
of lCen & Ostriken d 19991 ) and assumes that the temperature of IGM 
is 10 4 K at 2 > 3 and 10 6 /(1 + z) 3 ' 3 K at 2 < 3 and does not in¬ 
clude the contribution of massive clusters of galaxies. While these 
rough estimates cannot be compared quantitatively with the sim¬ 
ulation results, it demonstrates that these contributions are much 
smaller than those from hot gas in galaxy clusters in the simula¬ 
tion. 

The dashed gray line shows the analytical model of th e contri- 
butions of halos computed with the mass function of iTinker et al.1 
d20ld) and the Planck pressure profile with the mass bias of 1.2. 
The analytical model does a reasonable job describing the simula¬ 
tion result, although it is systematically lower than the simulation 
by 20% at 2 < 1. See section 4.2 for more detailed discussion on 
the analytical model. 

The dashed black lines show the scatter due to kSZ at 
150 GHz. The kSZ is most prominent at high redshift and starts 
to contribute even earlier than tSZ does. 

In figure [ 8 ] we show the evolution of the cumulative mean 
Compton Y. The solid lines show Y (> 2) = f j ’' 2 dz' dP/d 2 , 
while the dashed lines show Y(< 2 ) = J 0 “ dz' dY/dz. The latter 
clearly shows that the mean Compton Y receives significant con¬ 
tributions up to 2 ft; 2. We find that cluster-size halos with M > 
10 14 Mq /h contribute only 20% of the total signal, and nearly half 
of the signal, Y = 5 x 10“ 1 , co mes from M < 10 13 Mq /h (which 
can be detected by stacking; IPlanck Collaboration et alj 12013d ; 
iGralla et al.ll2014l : lGreco et al.ll2015l) and the diffuse barvons out¬ 
side halos. 


4 MEAN COMPTON Y 
4.1 Simulation results 


4.2 Analytical model 

The analy tical model for the m ean Compton Y is given analogously 
to Eq. [l2l ( lBarbosa et alJI 1996h 


The deep light-cone yields the mean Compton Y of Y = 1.18 x 
10 -6 (shown as the dotted horizontal line in figure^ for the cos¬ 
mological parameters used in the Magneticum Pathfinder simula¬ 
tion: 0 ‘s = 0.809 and f2 m = 0.272. We now study how Y builds 
up over cosmic time and how much objects of different masses con¬ 
tribute to Y. To this end, we proceed similar to the construction 
of the different slabs contributing to the light-cone. However, to 
avoid that the volume (and therefore the statistics) decreases with 


Y = 




Mm “ ,, dn(M,z) _ , 

M . dlnM ^i r w(M ’ zl 


(14) 


where yo is Eq.[l3]with l = 0. However, as the mean Compton Y 
receives significant contributions from lower mass halos compared 
to the power spectrum, we use M m in = 5 x 10 10 Mg/ft. As our 
goal in this section is to confirm the result of the previous section, 
we use the cosmological parameters of the Magneticum Pathfinder 
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z alog10(dY/dz*(z+1) 


Figure 8. Same as figure[7] but for the cumulative signals. The solid lines Figure 9. PDFs of dY/dln(l + z) at z = 0 (black), 0.5 (blue), 1 (red), and 

show Y(> z) = J.' 1 ' 2 dz' dY/d z' , while the dashed lines show Y(< 2 (green). 

z) = Jq 2 dz' dY /d z'. 


simulation. All the other details of the calculation are the same as 
in section [33.2l 

W e find Y = (0.78,0 . 83,0.99) x 1 0~ 6 for the mass func¬ 
tions of lBocquet et alj|20l5l) . ITinker et al.ld2008l) , and lTinker et all 
l20lol) . respectively. The mean Compton Y receives significant 
contributions from low-mass halos for which these mass functions 
differ significantly. In particular, the former two fitting functions 
do not satisfy the normalisation constraint on the mass function, 
/ 0 °° dM^oom dn/d In M200111 = p. A reasonable agreement be¬ 
tween the results from t he simulatio n and the analytical model with 
the mass function of iTinker et alj d20ld ). which does satisfy the 
normalisation constraint, is encouraging; however, more study on 
a low-mass end of the mass function is necessary. The remaining 
difference of order 20% relative to the simulation is due to the 
Planck pressure profile with the mass bias of 1.2 being slightly 
lower than the pressure profiles in the simulation in lower mass 
halos (M < 10 14 Mg). The same trend can be seen in the tSZ 
power spectrum at l > 3000 shown in figure [6] For example, we 
find? = (1.07,1.16) x 10 -6 with the mass biases of 1.15 and 1.1, 
respectively. However, the mass bias of 1.1 would yield too large a 
tSZ power spectrum at low multipoles, l < 1000, to agree with the 
simulation. 

Using the Planck 2015 parameters, the mass function of 
ITinker et all d2Q10h and the mass bias of 1.2, we find Y = 
1.32 x 1CU 6 . Using this to scale the simulation result, we find 
? = 1.57 x 1 0~ 6 . _ 

Recently, |HjU_etal3 l|2015| ) use t he analytical model with the 
mass function of Tinker et aiT hoiol) and the pressure profile of 
iBattaelia et all d2012h to obtain ? = 1.58 x 10 -6 f or the WMAP9 
param eters with erg = 0.817 and fl m = 0.282 dHinshaw et al.1 
120131) . A few minor details on the calculation are different: while 
they integrate the pressure profile out to 2.5 times the virial ra¬ 
dius, we integrate out to 6rsoo; their redshift integration is from 
2 = 0.005 to 6; and the lower boundary of their mass integra¬ 


tion is M m in = 5 x 10 5 Mq/Zi. Changing the upper integration 
boundary from a; max = 6 to 2.5r v i r /rsoo in Eq.[l3]and adjusting 
the integration boundaries in Eq. |T4] we find ? = 1.11 x 1(P 6 
for the Planck pressure profile with the mass bias of 1.2 and the 
same WMAP9 parameters. This value is significantly lower than 
their value. We find that this is due to the difference in the pressure 
profiles; while the pressure profiles of the Magneticum Pathfinder 
simulation, the Pl anck pressu re profile with the mass bias of 1.2, 
and the profile of lBattaglia et al.l 1120121) agree i n the high-mass end, 
M 500 > 2 x 10 14 Ma\/h Jpianck Collaboration et al .lllO 1 3~3) , the 
latter profile gives significantly larger pressure than the Planck pro¬ 
file in lower masses that dominate in Y. 

These discrepancies need to be understood before we obtain 
an accurate estimate of the expected level of Y. Nevertheless, the 
original conclusion from the previous generation o f cosmological 


hydrodynamical simulati ons |Rgfregigr_et_al. 2000; d a Silva et al.l 


120001 ; ISeliak et~akl 1200 ll ; IZhang etak 2004b ) seems robust: the 
expected Y from the large-scale structure is of order 10 -6 and 
is only one order of magnitude lower than the FIRAS bound. It 
is also encouraging that none of these estimates exceed the new 
Planck bound on the fluctuating pa rt of the mean Y parameter, 
? < 2.2 x 10 -6 khatn & Sunvaevll2015l) . 


4.3 Buildup of the Compton Y PDF 

Finally, we study contributions from different redshifts to the build 
up of the PDF of the Compton Y signal. Figure[9]shows the PDF of 
dY/dln(l + z) at four different redshifts. Although these PDFs look 
at first glance similar to the one of the full light-cone (as shown in 
figure |5J, the tail does not follow a simple power law. In general, 
with decreasing redshift, the PDF becomes broader and less sharply 
peaked, especially when compared to the PDF at 2 = 2. 

In figure [TO] we show PDFs with high-mass halos above a 
certain mass threshold removed from the simulation. At any given 
time, the tail is dominated by the most massive halos of that time. 
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2=0.0 z = 0.5 




z= 1.0 2 = 2.0 




Figure 10. PDFs of dY/dln(l + z) at 2: = 0 (top left), 0.5 (top right), 1 (bottom left), and 2 (bottom right). In each panel, we show PDFs with high-mass 
halos above a certain mass threshold removed from the simulation. The colours indicate the same mass thresholds as shown in figure[7] 


5 CONCLUSIONS 


The Magneticum Pathfinder simulations, state-of-the-art cosmo¬ 
logical and hydrodynamical simulations, follow in detail the ther¬ 
mal and chemical evolution of the ICM as well as the evolution 
of super-massive black holes and their associated feedback pro¬ 
cesses. These simulations reproduce the avera ge ICM p ressu re pro¬ 
files measured by Planck dPlanck Collaboration et all 120135) and 


SPT dMcDonald et al.ll2014l) . At the same time, the stellar mass 
functions of galaxies and the luminosit y functions of the AGN 
population agree well with observations dHirschmann et al.ll2014l ). 
The improved numerical methods and increased computing power 


available today have enabled these simulations to follow a large 
enough cosmological volume to construct realistic light-cone maps, 
which can be compared with observations in detail. 

In this paper, we have computed the tSZ and kSZ effects to¬ 
ward the counterpart of the Coma cluster in the local universe sim¬ 
ulation, and statistics of the SZ effects from the light-cone maps, 
including the one-point PDF and power spectrum of tSZ and kSZ, 
and the mean Compton Y parameter. We have then compared these 
predictions on tSZ with the Planck, SPT, and ACT data. Our find¬ 
ings are summarised as follows: 

• The tSZ radial profile of Coma in the local universe simulation 
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embedded in the background from the deep light-cone agrees well 
with that in the Planck data. 

• The local universe simulation predicts that the halo of Coma 
is moving away from us at « 400 km/s with respect to the CMB 
rest frame, thus yielding a negative kSZ within the central region 
of Coma. The magnitude of kSZ is ten percent of tSZ at 150 GHz. 
On the other hand, because Coma is a merging system, we find a 
significant relative motion of the core (even increasing the negative 
signal in the center) and a significant positive kSZ in the outskirt, 
which comes from a infalling sub-structure moving toward us. This 
makes a positive kSZ contribution to tSZ at 150 GHz at distances 
beyond 1 Mpc from the center of Coma. 

• The predicted one-point PDF of the Compton Y agrees with 
that measured by Planck, once the simulations are smoothed to the 
resolution of Planck’s Y map (10 arcmin). Given the much smaller 
beam size, we expect that ACT- and SPT-like instruments will see 
almost the full PDF that is resolved by our simulations. The tail of 
the full PDF follows a power-law with an index of —3.2. 

• The tSZ power spectrum measured from the simulation 
agrees with that of the Planck data at all multipoles up to l ~ 
1000, once the power spectrum is rescaled to Planck 2015’s 
“TT+lowP+lensing” cosmological p aram eters with — 0.308 
and as = 0.8149 dPlanck Collaboration et alj|2015al) . We have 
confirmed and understood this result using the analytical model. 

• Consistent with the previous work, we continue to find the pre¬ 
dicted tSZ power spectrum at l = 3000 that is significantly higher 
than that estimated by ACT and SPT. Whether this poses a chal¬ 
lenge to theory is unclear, but our prediction is still well below the 
firm upper bound on tSZ given by the SPT data points with the 
primary CMB subtracted. 

• The simulation predicts the mean Compton Y value of 1.18 x 
10 -6 for f l m = 0.272 and as = 0.809. When the contributions 
from halos above a virial mass of 10 13 M@/h are removed, we 
find Y = 5 X 10 - '; thus, nearly half of the signal comes from such 
low-mass halos and diffuse gas outside halos. This remaining signal 
would pose a challenge to detecting the primordial y-distortions. 

• Using the analytical model, we scale the Compton Y value 
from the simulation to the Planck 2015 parameters with = 
0.308 and cr 8 = 0.8149, finding Y = 1.57 x 10 -6 . This is still 
lower than, but not far away from, t he new Planck bound, Y < 
2.2 x 10~ 6 jKhatri & Sunyaev! 20151) . 

• The one-point PDF and the power spectrum of kSZ from our 
simulations agree broadly with the previous work. While our box 
size is large, the contribution to kSZ is still dominated by the largest 
modes within the box and originates mainly from high redshifts. 
Therefore, unlike for tSZ, we have not yet obtained a reliable, con¬ 
verged result on kSZ on large scales, l < 1000. Simulations fol¬ 
lowing even large cosmological volumes are needed. 


In short, the main conclusion from our study is that all the 
properties of tSZ found in the Magneticum Pathfinder simula¬ 
tion and the local universe simulation agree well with the Planck 
data. This includes the tSZ power spectrum, which was previ- 
ously found to be in tension with the Planck 2013 par ameters 
dPlanck Collaborati on et alj l2014bl : Mc Carthy et al.l 120141 ). Now, 
the tSZ power spectrum calculated for the Planck 2015 parameters 
including CMB lensing information agrees with the measurement 
at all multipoles up to l fa 1000. 
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